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Errata  to: 

MRC  Technical  Summary  Report  #2915 

Zero  and  Negative  Masses  In  Finite  Element  Vibration  and  Transient  Analysis 
David  S.  Malkus  and  Michael  E.  Plesha 


page  1,  paragraph  2 
page  9,  paragraph  1 
page  16,  equation  (37) 

page  16,  equation  (40) 
page  18,  paragraph  5 

page  21,  paragraph  2 
page  22 


negtlve  should  be  negative 

0  £  1  <_  ...  should  be  1  £  1  _<  . . . 

last  term  In  equation  should  be  postmultlplled 
by  Instead  of  {v  } 

4/At^  should  be  At^/4 

last  sentence:  reference  to  equation  (46)  should 
be  to  equation  (47) 

an  viable  should  be  a  viable 

somes  sense  should  be  some  sense 


Section  6  following  references  on  original  work  in 

partitioned  transient  analysis  should  be  added: 

T.  Belytschko  and  R.  Mullen,  Mesh  Partitions  of 
Explicit-Implicit  Time  Integration,  in  Formulations 
and  Computational  Algorithms  in  Finite  Element 
Analysis,  K.  J.  Bathe,  et  al  eds.,  MIT  Press,  1976. 

T.  Belytschko  and  R.  Mullen,  Stability  of 
Explicit-Implicit  Mesh  Partitions  in  Time  Integration 
Int'l.  J.  Numerical  Methods  in  Engr.,  12, 

1575-1586  (1978). 
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Computer  in  Applied  Mechanics  and  Engineering. 


UNIVERSITY  OF  WISCONSIN-MADISON 
MATHEMATICS  RESEARCH  CENTER 


ZERO  AND  NEGATIVE  MASSES  IN  FINITE  ELEMENT  VIBRATION 
AND  TRANSIENT  ANALYSIS 
David  S.  Malkus*  and  Michael  E.  Plesha** 

Technical  Summary  Report  #2915 
February  1986 

ABSTRACT 

^  Mass  matrix  lumping  by  quadrature  is  considered.  Accuracy  requirements  seem  to 
dictate  the  use  of  zero  or  negative  masses  for  multi-dimensional  higher-order  elements. 
It  is  shown  that  the  zero  and  or  negative  masses  do  not  destroy  the  essential  algebraic 
properties  of  the  discrete  eigenproblem,  in  spite  of  the  negative  or  infinite  eigenvalues  which 
may  result.  Explicit  transient  methods  require  positive  definite  lumping  which,  for  some 
elements,  may  only  be  achieved  by  sacrificing  accuracy  to  avoid  the  negative  or  zero  masses 
that  would  render  the  lumping  indefinite.  An  implicit-explicit  time  integration  method 
based  on  quadratic  triangles  with  optimal  lumping  is  devised,  analysed,  and  tested.  It 
treats  the  nodes  with  nonzero  masses  explicitly  and  the  nodes  with  zero  masses  implicitly. 
Analysis  and  numerical  tests  show  that  this  formulation  is  optimally  accurate  and  less 
costly  than  a  similar  method  with  nonzero  masses,  based  on  optimally  lumped  biquadratic 
rectangles.  The  method  is  also  found  to  be  substantially  more  accurate  than  the  fully 
explicit  method  based  on  lumping  the  triangular  elements  in  an  ad-hoc  fashion  to  retain 
non-zero  masses.  /  • 
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ZERO  AND  NEGATIVE  MASSES  IN  FINITE  ELEMENT  VIBRATION 
AND  TRANSIENT  ANALYSIS 

David  S.  Malkus*  and  Michael  E.  Plesha** 

1.  INTRODUCTION 

Mass  matrix  lumping  is  a  technique  in  finite  element  transient  and  vibration  analysis 
whereby  the  banded,  symmetric,  positive  definite  matrix  representing  the  inner  product 
(i.e.,  the  consistent  mass  matrix)  is  replaced  by  a  diagonal  matrix  which  is  in  some  sense 
equivalent.  This  is  done  for  two  basic  purposes:  first,  to  make  truly  explicit  finite  element 
transient  analysis  possible,  and  second  to  save  computer  storage  space  and  operations  in 
the  solution  of  eigenproblems  associated  with  the  finite  element  vibration  problem.  There 
are  two  additional  benefits  associated  with  lumped  masses:  in  transient  analysis,  lumping 
lowers  the  highest  wave  speed  in  the  mesh  significantly,  with  an  attendant  raising  of  the 
critical  time  step  for  explicit  integration  schemes;  in  eigenanalysis.  lumping  removes  the 
strict  upper-bounding  of  the  discrete  eigenvalues.  The  result  can  be  that  the  eigenvalues 
remain  optimally  accurate  in  terms  of  convergence  rate,  but  have  a  better  constant  in  the 
error  bound. 

The  role  of  lumping  in  transient  analysis  makes  it  an  essential  tool  in  finite  element 
analysis.  Many  problems  in  structural  analysis  are  nonlinear  hyperbolic  problems  which 
can  be  treated  explicitly.  The  geometric  complexity  of  many  practical  problems,  such  as 
in  reactor  technology  requires  so  many  finite  elements  that  any  fully  implicit  treatment 
of  the  nonlinearity  would  be  impractical.  Thus  the  basic  idea  behind  lumping  stems  from 
practical  necessity  and  has  a  ready  physical  interpretation:  the  continuous  distribution  of 
mass  in  the  body  is  replaced  by  one  in  which  the  mass  is  concentrated  at  the  nodal  points. 
That  this  has  the  desired  effect  on  the  mass  matrix  can  easily  be  seen  by  replacing  the 
mass  density  function  by  the  sum  of  delta  functions  with  support  at  the  nodes.  For  low- 
order  elements,  the  requirements  of  preserving  the  total  mass  of  the  element  and  respecting 
the  symmetry  of  the  distribution  about  natural  axes  of  the  element  uniquely  determines 
the  lumped  distribution.  For  higher-order  elements,  many  distributions  seem  to  respect 
natural  symmetry,  but  most  seem  to  lead  to  schemes  with  reduced  accuracy. 

An  approach  to  recovering  optimal  accuracy  with  lumping  schemes  is  to  require  each 
element  to  exactly  reproduce  as  many  moments  of  the  mass  distribution  as  possible.  This 
turns  out  to  be  equivalent  to  the  use  of  a  nodally  based  quadrature  formula  to  evaluate 
the  L2  inner  product.  This  is  the  fundamental  observation  of  ref.  1.  Once  we  recognize 
that  lumping  and  numerical  quadrature  correspond,  then  we  can  be  guided  by  the  error 
analysis  of  the  latter  in  devising  lumping  schemes.  Here  we  consider  only  second-order 
problems,  and  in  that  case,  ref.  1  argues  that  we  should  choose  a  nodal  formula  of  degree 
of  precision  2p  -  2  where  p  is  the  degree  of  the  maximal  complete  polynomial  on  each 
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element.  When  that  is  done,  the  error  introduced  by  quadrature  will  be  bounded  by 
C/i^  l|u  and  thus  be  on  the  same  order  as  the  overall  error  of  the  scheme,  and  no 
accuracy  will  be  lost. 

One  quickly  finds  that  the  degree  of  precision  requirement  for  retaining  optimal  ac¬ 
curacy  is  not  a  trivial  one.  In  fact,  for  elements  with  p  >  3.  the  traditional  equal  spacing 
of  nodes  in  the  reference  element  will  no  longer  suffice  jlj.  A  multidimensional  generaliza¬ 
tion  of  Lobatto  2j  integration  is  called  for;  vertex  nodes  must  be  included,  and  boundary 
nodes  must  remain  on  boundaries,  but  the  nodal  locations  must  be  chosen  so  as  to  oth¬ 
erwise  maximize  the  degree  of  precision  For  one-dimensional  problems,  the  Lobatto 
weights  can  be  shown  to  be  positive  [2,.  This  is  important,  since  in  all  dimensions  the 
quadrature  weights  are  the  diagonal  entries  of  the  lumped  mass  matrix.  But  in  two  or 
more  dimensions,  a  fundamental  difficulty  in  numerical  analysis  arisies:  maintaining  posi¬ 
tive  weighting  and  specifying  quadrature  point  locations  can  be  mutually  exclusive.  Two- 
dimensional  elements  for  which  this  conflict  occurs  are  schematically  illustrated  in  Figure 

1.  The  appearance  of  zero  and  negtive  weights  seems  to  leave  the  practitioner  four  choices: 
(1)  use  low  order  elements,  (2)  use  elements  which  are  tensor  products  of  one  dimensional 
elements,  (3)  use  lumping  schemes  which  sacrifice  degree  of  precision  to  the  maintenance 
of  positive  weights,  or  (4)  employ  a  non-positive-definite  mass  matrix  with  coefficients  that 
can  be  positive,  zero  or  negative.  The  purpose  of  this  paper  is  to  demonstrate  that  (4)  is 
a  viable  computational  alternative. 

2.  THE  NONPOSITIVE  MASS  PROBLEM 

There  appear  to  be  at  least  three  difficulties  with  nonpositive  mass  lumpings,  beyond 
the  obvious  physical  peculiarity.  Actually,  the  physical  peculiarity  cannot  be  a  particularly 
serious  problem,  because  the  lumping  scheme  is  devised  to  produce  optimally  accurate 
eigenvalues.  Their  presence  in  the  spectrum  and  their  accuracy  can  thus  be  guaranteed 
a-priori,  and  any  other  modes  present  in  the  problem  must  take  their  place  within  the 
unphysical  portion  of  the  spectrum,  which  is  there  whether  or  not  ^he  masses  are  lumped. 
But  there  are  a  number  of  possible  technical  difficulties.  In  eigenanalysis,  zero  masses  are 
known  not  to  be  problematical,  and  can  in  fact  be  put  to  good  use  in  static  condensation 
[3).  But  in  transient  analysis  zero  masses  seem  to  defeat  the  very  purpose  of  lumping.  An 
explicit  multistep  method  would  have  the  form 

m{Dn}  =  {T({Dn-l}ADn-2} . {Dn- k} ,  {Fn-l}))  (l) 

Clearly  no  unique  solution  can  exist  when  \M]  is  singular,  even  if  diagonal.  Negative 
masses  seem  even  more  troublesome.  The  finite  element  eigenproblem 

\K\{D}  =  X{M]{D}  (2) 

could  have  negative  eigenvalues  when  |Afj  is  indefinite.  For  statically  well-determined 
problems  we  shall  presume  (A’j  is  positive  definite.  In  that  case  it  is  clear  that  there  will 
be  as  many  negative  eigenvalues  as  negative  masses.  However,  for  free  body  vibrations, 
'.K\  is  singular.  In  such  a  case  there  may  or  may  not  be  negative  eigenvalues  when  \M' 
has  negative  entries.  One  cannot  be  sure  a-priori  because  when  M]  is  indefinite  and  j/f  j 
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FIGURE  1 

Optimal  lumpings  for  some  common  1-  and  2-D  elements.  Tri¬ 
angular  elements  are  equilateral  reference  elements  in  area  co¬ 
ordinates.  Complete  symmetry  of  lumping  formula  in  area 
co-ordinates  is  assumed,  determining  nodal  locations  and  un¬ 
specified  weights  completely  from  given  information.  Square 
elements  are  isoparametric  master  elements.  Complete  sym¬ 
metry  in  Cartesian  co-ordinates  uniquely  determines  lumping 
formulas.  “A”  denotes  area  of  2-D  elements. 
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is  singular,  the  fundamental  structure  of  the  eigensystem  can  be  quite  different  from  the 
case  in  which  either  matrix  is  positive  definite.  A  simple  example: 


\K]  = 


\M] 


-1  0 

0  1 


Let 


Then  the  matrix  pencil  \K'<  -  A'Af ;  is  equivalent  to  4 


-X  1 
0  -A 


(3) 


We  therefore  have  a  multiple  zero  eigenvalue  and  a  nonlinear  divisor.  No  negative  eigenval¬ 
ues  exist,  auid  the  matrix  pencil  is  no  longer  diagonalizable.  Examples  can  be  constructed 
in  which  the  eigenvalues  are  complex.  Even  worse,  if  there  are  zero  masses,  examples  can 
be  found  in  which  there  fails  to  be  a  complete  set  of  eigenvalues  [4].  The  simplest  way  this 
can  happen  is  for  {K\  and  iM]  to  have  a  common  null-space. 

The  above  example  is  clearly  contrived.  It  would  not  arise  in  finite  element  analysis; 
[K]  is  a  simple  stiffness  matrix  for  a  linear  element,  but  the  mass  matrix  is  lumped  in  a 
way  that  does  not  preserve  mass.  To  put  it  another  way,  the  quadrature  formula  whose 
weights  are  the  diagonal  entries  would  violate  the  degree  of  precision  requirement:  the 
degree  of  precision  is  not  even  zero,  whereas  2p  -  2  =  0.  The  fear  in  using  negative  masses 
is  that,  while  quadrature  accuracy  may  force  the  existence  of  accurate  eigenvalues,  they 
may  be  tied  up  in  nonlinear  divisors,  have  the  wrong  multiplicity,  and  be  hard  to  compute. 
There  could  conceivably  be  negative  eigenvalues  of  the  same  magnitude  of  the  accurate 
eigenvalues  leading  to  spurious,  unstable,  low  frequency  modes. 

The  first  task  to  which  we  turn  our  attention  is  to  generalize  what  we  observed  in 
the  simple  example,  the  fact  that  preservation  of  mass  and  higher  moments  of  mass  rules 
out  the  pathologies.  But  this  is  easier  to  tackle  in  an  equivalent  form;  satisfaction  of  the 
degree  of  precision  requirement  rules  out  all  the  difficulties  allude  to  above. 


3.  MASS  MATRIX  LUMPING  BY  QUADRATURE 

We  shall  deal  with  second-order,  self-adjoint  problems  in  which  the  natural  modes  of 
the  exact  problem  make  the  Rayleigh  quotient,  Q(u),  stationary, 


Q(u)  = 


a(u,u) 

6(u,u) 


(4) 


where  a(-,  •)  represents  the  strain  energy,  6(-,  •)  represents  the  L2  inner  product,  and  for 
some  real  '1 


a(u,u)  > 

for  u  e  ^  =  H.  denotes  the  space  of  functions  with  square  integrable 

first  partials  satisfying  essential  boundary  conditions  appropriate  to  the  vector  compo¬ 
nent  of  the  solution.  For  plane  or  1  -  D  problems,  from  which  our  examples  are  drawn,  we 


may  consider  that  =  0  for  one  or  two  factors  of  the  trial  space.  Generalization  of  our 
results  and  examples  to  3  •  D  is  straight-forward. 

It  will  be  convenient  here  to  use  fairly  standard  engineering  finite  element  notation 
[3i;  u  €  H  will  be  approximated  in  element  “e”  by 

=  (5) 


where  {d}  represents  the  nodal-value  degrees  of  freedom  and  \Ne]  is  the  element  shape- 
function  matrix  which  has  the  form 


Nu  0  ...  Sme  0 

0  N  If  ...  0  A^me 


(6) 


for  an  m-noded  element  in  2  -  D.  Generalizations  to  1  -  and  3-D  problems  are  c  bvious. 
Note  that  equation  (6)  implies  a  nodal  ordering  in  which  the  degrees  of  freedom  associated 
with  a  given  node  are  grouped  together  in  {d}  (3).  Element  shape  functions  are  pieced 
together  in  the  usual  manner  to  produce  a  global  nodal  representation 

u**  =  [N]  {D}  (7) 


where 


Ni 
0  Ni 


Nm  0 
0  Nm 


(8) 


which  illustrates  the  M-noded,  2  -  D  case.  Satisfaction  of  the  essential  boundary  conditions 
in  equation  (7)  is  assumed  to  be  accomplished  by  specification  of  appropriate  entries  of 
{£>}  to  prescribed  values  [3]. 

The  free-structure  stiffness  matrix  is  assembled  from  element  matrices 


m  =  jjB.f  \E]{B,]  dV 


(9) 


where  [E\  is  a  symmetric  matrix  of  material  coefficients  and 


(10) 


[Pi  is  a  matrix  of  first-order  partial  differential  operators  relating  stress  to  strain  (or  strain- 
rate  in  fluid  problems).  For  example,  in  plane  elasticity 


1P1  = 


0 

0 

A 

>?y 

A 

dz 

(11) 


Boundary  conditions  implied  by  are  assumed  to  be  homogeneous  and  are  enforced 
by  omission  of  corresponding  rows  and  columns  during  assembly  l3j.  Entries  in  the  global 
nodal-value  vector.  {D),  are  omitted  along  with  corresponding  columns  in  |A^j.  The  space 


of  all  such  functions  defined  by  equation  (7)  is  denoted  by  S*.  The  global  stiffness  matrix, 
possibly  dimensionally  reduced  by  imposing  boundary  conditions,  is 


where 


'K]  =  [  EHB]  dV 
hi 


B]  =  [PMJVi 


(12) 


(13) 


For  future  reference  in  dynamic  problems  we  note  that 


{<,}  =  |B,|  {d) 


and 


(<)  =  |B|  {D} 


(15) 


Equations  (14)  and  (15)  give  the  virtual  work  done  against  internal  stresses  by  the  displace¬ 
ments.  The  matrix  of  a(',  •)  in  equation  (4)  is  [/f  j.  With  no  essential  boundary  conditions, 
"](  =  0. 

The  mass  matrix  gives  the  matrix  of  6(',  •)  in  equation  (4)  and  is  defined  on  the  element 
and  global  level  by 


\m,]  =  j^p\N,]^{N,]  dV 
[M]  =  f  dV 


(16) 


where  p  is  the  material’s  mass-density  function.  In  what  follows,  p  may  be  non-constant, 
but  p  >  0  is  assumed.  When  the  integrations  in  equation  (16)  are  exact,  [me]  and  \M\ 
are  called  “consistent”.  To  lump  the  mass  matrix  by  quadrature,  the  nodes  are  located  in 
order  to  maximize  the  degree  of  precision  of  a  nodal  quadrature  formula,  while  retaining 
interelement  compatibility  ill,  and  equation  (16)  is  integrated  with  the  nodal  quadrature 


rule,  resulting  in 


rrie 


0  0] 

0  0 


m 

1=1 


1  0 
0  1 


0  0...  ]  0  ...  0  0 

0  0  ...  0  ]  ...0  0 

L 


0  0 
0  0 


PiWiJi  j72j 

P2W2J2  (/2i 


o 


(17) 


O  • 

L  PmVJmJm\h\^ 

where,  again,  an  m-noded  2  -  D  element  illustrates  the  generality;  Pi  and  J,  are  evaluations 
of  the  density  and  any  isoparametric  or  curvilinear  Jacobian  implied  in  “dK”  at  the  element 
nodes,  respectively,  and  the  are  the  quadrature  weights.  1/2]  is  the  2x2  identity.  In  this 
notation  superposed  tilde  denotes  quantities  evaluated  by  numerical  quadrature.  Clearly 
the  [me]  can  be  assembled  into  a  global,  lumped  mass  matrix,  [M],  which  is  diagonal. 


4.  LUMPING  AND  EIGENSTRUCTURE 

The  basic  result  we  use  to  investigate  the  effect  of  lumping  on  the  pencil  \K]  —  \\M] 
is  that  of  G.  Fix  {5j,  which  deals  with  the  accuracy  of  applying  a  quadrature  formula  to 
equations  (9)  and  (16)  (on  the  element  level)  to  produce  numerically  integrated  forms  c(-,  ■) 
and  6(-,-)  when  numerically  integrated  stiffness  and  mass  matrices  are  assembled.  The 
analysis  applies  to  lumping  by  quadrature.  The  major  result  is  bsised  on  the  satisfaction 
of  two  conditions; 

(A)  Stability;  Let  u**,  v**  € 


sup 


c(u*’,v*')  -  c(u'*,v 


h 


h  11  =1 


(B)  Accuracy;  Let  u**  G  5^  be  a  finite  element  eigenfunction  and  v**  G 


(18) 


|j  jj  ^  j  I  -  (19) 

where  c(-,-)  denotes  o(-,-)  or  t>(-,-),  and  h  is  the  mesh  parameter,  typically  the  maximal 
element  diameter. 

Let  u**  denote  eigenfunctions  which  are  stationary  points  of 


Note  that  a  =  a  is  possible  (mass  matrix  quadrature  only),  or  a(-,-)  could  be  subjected 
to  quadrature  as  well  (the  same  or  different  from  that  applied  to  6(-,  •)).  Let  A/,  and  \h 
denote  eigenvalues  of  Q  and  Q,  respectively.  Fix’s  major  result  states  that  if  conditions 
(.4)  and  (B)  are  met  for  a  given  u**,  then  there  is  a  u**  such  that 

u**  -  .  j  <  -  C'Sl!  =  6  (21) 

for  constants  C  and  C,"  independent  of  h.  For  the  symmetric  problems  considered  here 
a  simple  argument  leads  to  the  conclusion  that  A/,  -  A;,:  =  O(S^)  for  A;,  =  <?(u*‘)  and 
Xh  -  Qlu**).  A  simple  consequence  of  condition  (A)  is  that  if  it  is  satisfied  and  7  >  0, 
there  is  an  h.  such  that  for  all  l\  <  h,,.  a(v**,v**)  >  ,  v**  ^  with  ')  >  0.  This  will  soon 

play  a  crucial  role  here. 

In  order  to  make  -  0(/?'’  I  u  ^  j)  where  p  is  the  degree  of  the  maximal  complete 
polynomial  on  each  element  and  u  is  a  stationary  point  of  equation  (4),  the  degree  of 
precision  of  the  quadrature  formula  must  be  sufficiently  high.  If  it  is,  then  no  inaccuracy 
will  be  incurred  by  the  quadrature  which  is  of  lower  order  in  h  than  ||  u  -  u**  ||^.  In  ref. 
1  the  statement  is  made  that  the  degree  of  precision  should  be  2p  -  2.  According  to  the 
analysis  in  ref.  5,  this  is  both  necessary  and  sufficient  for  non-isoparametric  triangles, 
which  will  be  the  elements  of  most  concern  to  us  here.  The  degree  of  precision  criterion 
given  in  jS]  is  more  stringent  than  2p  -  2  for  rectangles  and  appears  to  be  excessively  so, 
to  the  point  of  excluding  optimal  lumping  formulas  based  on  product  Lobatto  formulas 
for  bi  -  p  degree  rectangles  (see  Figure  1).  The  reason  is  that  Fix’s  analysis  applies  to 
integration  of  the  stiffness  and  mass  matrices,  whereas  if  only  the  mass  matrix  is  perturbed. 
Fix’s  result  can  be  sharpened.  In  the  Appendix  to  this  paper  we  provide  a  proof  which 
applies  to  undistorted  bilinear  and  biquadratic  elements  and  shows  that  (A)  and  (B)  are 
satisfied,  whereas  the  degree  of  precision  is  higher  than  2p-  2  and  lower  than  required  by 
Fix. 

In  general  the  authors  have  found  that  optimizing  the  degree  of  precision  subject  to 
the  appropriate  nodal  constraints  leads  to  stable  and  optimally  accurate  schemes.  We 
know  of  no  completely  general  proof  of  this  which  applies  for  non-trivial  isoparametric 
transformations.  But  computational  experience  seems  to  suggest  that  optimal  lumping  by 
quadrature  produces  accurate  schemes  for  isoparametric  elements.  For  the  elements  we 
use  here  our  own  proof  applies.  For  those  elements  for  which  no  optimality  proof  is  known, 
it  seems  to  stand  to  reason  that  choosing  the  lumping  formula  with  highest  possible  degree 
of  precision  is  the  wisest  practical  choice  from  the  point  of  view  of  retaining  accuracy. 

We  turn  our  attention  to  the  eigenstructure  of  the  lumped  pencil,  [/f]  —  A[A/].  We 
first  consider  the  constrained  body  case  in  which  |/C)  is  positive  definite.  From  ref.  4  it 
follows  immediately  that  there  is  a  nonsingular,  real  [Ql  for  which 


-  A  Mi):Q' 


■  Ql  -  A/?i 

O2  —  A^2 


o 


o 


(22) 
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Clearly  a,,/?,  are  real  and  a,  >  0.  Assume  the  ordering  is  such  that  the  k  zero  3i  are 
ordered  first,  so  /3i  =  0,  0  <  i  <  k  <  n.  A  few  simple  facts  should  be  apparent: 

(1)  The  eigenvalues  are  aif3i,  k  <  i  <  n  and  oo  {k  times). 

(2)  If  Ao  is  any  real  number,  the  eigenvalues  of  \K'  -  (A  -  Ao)[Adr]  are  a,/j3,  +  Ar,.  k  <  i  <  n 
and  oc  {k  times). 

We  now  turn  our  attention  to  the  case  in  which  \K  has  rigid  body  modes  and  is  only 
semi-definite. 

Lemma:  Let  Condition  (A)  hold  for  the  lumping  formula.  Given  any  Ao  >  0,  there  is  an 
h  sufficiently  small  such  that  for  all  h  <  h,,.  K ]  -  A,,  A/  is  positive  definite. 

Proof:  \K]  +  XfJM]  is  a  positive  definite  stiffness  matrix.  A'  -rAci[M]  is  a  stable  quadrature 
evaluation  of  the  consistent  matrix.  Stability  implies  the  existence  of  a  7  such  that 

c(u,u)  -r  Ao6(u.u)  >  7  ii  U  \'‘fj 

I 

Theorem  1:  For  sufficiently  small  h  there  exists  \Q  sis  in  equation  (22)  diagonalizing  the 
pencil  [A”]  -  A[Mj  whether  \K\  is  positive  definite  or  semi-definite,  whenever  the  lumping 
quadrature  is  stable. 

Proof:  Pick  any  Ao  >  0.  Let  h  be  sufficiently  small  so  that  [A]  +  Ao[M]  is  positive  definite. 
Equation  (22)  applies  to  ([A]  -h  Ao[M))  -  A[M]  =  [A]  -  (A  -  Ao)[M).  As  observed  above, 
reshifting  by  -Ao[M]  does  not  change  the  eigenstructure  and  results  in  [A]  -  A[A^j.  | 
Remarks: 

(1)  Since  Ao  can  be  arbitrarily  small  in  Theorem  1,  it  is  no  surprise  that  in  practice  the 
result  is  observed  to  hold  without  making  h  small  amd  even  to  hold  for  one  element. 

(2)  Theorem  1  says  that  the  divisors  of  the  pencil  (A)  -  A[a5^|  are  linear  and  furthermore 
the  pencil  is  diagonalizable  by  |M]-orthogonai  transformations. 

(3)  [Q]‘^[M\  [Q\  has  the  same  signature  as  [M],  so  that  the  number,  fc,  of  zero  3i  is  the  same 
as  the  number  of  zero  masses.  The  number  of  negative  0i  (and  thus  negative  eigenvalues) 
is  the  same  as  the  number  of  negative  masses  (assuming,  of  course,  h  <  ho). 

(4)  There  is  no  restriction  on  the  number  of  zero  and  negative  masses,  as  long  as  the 
stability  condition  holds.  Our  stability  proof  requires  degree  of  precision  is  at  least  zero 
(see  Appendix).  If  we  let  {1}^  represent  a  vector  with  all  entries  associated  with  the  first 
degree  of  freedom  at  each  node  equal  to  1: 

>  Pminvol.  (n,)  ^  0 

{1}^(M1{1}  >p,ni„vo/.(n)  #0 

where  17,  =  element  c,  and  f)  =  problem  domain,  and  Pmin  is  the  minimum  value  of  p  in 
the  appropriate  domain.  Thus  the  pathologies  of  the  example  of  the  previous  section  are 
ruled  out.  and  elements  wiht  nonpositive  mass  cannot  result  from  stable  lumpings. 

Next  we  consider  the  negative  eigenvalues.  The  worry  is  that  they  could  lead  to 
unstable  modes  with  low  enough  frequency  to  represent  mode  shapes  with  significant  ap¬ 
proximation  value.  This  cannot  happen,  at  least  on  a  fine  enough  mesh: 

Theorem  2:  For  a  stable  lumping,  given  any  N  >  0,  there  is  a  sufficiently  small  ho  =  ho(N) 
such  that  for  all  h  <  ho  any  a,;  d,  <  0  satisfies  |q,//3,|  >  N. 


'4u4^ibJ 


Proof:  Choose  ho  so  that  for  all  h  <  ho,  (/f|  + A'[A/j  is  positive  definite.  Then  a,  +N0i  >  0, 
and  if  at//?t  <  0  this  must  happen  because  0i  <  0{oLi  >  0).  Thus 

a.7/3.  <  -N 


-a,  13,  '>  N 

I 

Corollary  2.1:  Under  the  conditions  of  Theorem  2  the  pencil  [A'j  -  A[Mj  is  always  regular; 
that  is  there  are  no  eigenvalues  of  the  form  “0/0”. 

Proof:  In  the  above  proof  lake  N  >  0  arbitrarily  small,  a,  +  3iN  >  0.  This  is  a 
contradiction  if  a,  =  3,  -  0.  I 

Corollary  2.2:  Under  the  conditions  of  Theorem  2,  there  is  a  complete  set  of  eigenvalues 
of  |K]  -  AjM],  namely,  at,  3,.k  <  i  <  n  and  oo  {k  times). 

Remarks: 

(1)  Experiments  show  that  the  negative  eigenvalues  are  large  in  magnitude  relative  to  the 
lower  spectrum,  even  when  h  is  not  small. 

(2)  The  unstable  modes  due  to  a  negative- weighted  lumping  formula  are  always  mesh- 
dependent  and  are  driven  higher  and  higher  in  frequency  with  mesh  refinement. 

(3)  The  results  of  the  Theorem  and  Corollaries  flow  from  stability  and  the  fact  that  for 

small  enough  h,  \K\  +  can  be  made  positive  definite. 

This  may  seem  hard  to  believe  in  that  7V|A/]  adds  negative  weights  to  the  diagonal  of 
j/Cj.  But  closer  inspection  and  ref.  |6]  shows  that  the  diagonal  of  |/f]  +  N\M]  has  entries 
with  the  following  orders  of  h: 

1  -  D:  \K]  =  0{l/h),{M]  =  0{h).  Thus  diag{\K\  -  N|M])  -  0{l/h)  ±  0{Nh) 

2-  D:  \K]  =  0(1),  [M]  =  0{h^).  Thus  diag{{K]  -I-  N\M])  -  0(1)  ±  0{Nh^) 

3- D:  {K\  =  0(h),  [M\  =  0(h%  Thus  diaff({/fl  +  A^IM])  -  o\h)  ±  0{Nh^). 

In  the  light  of  the  above,  the  result  is  less  surprising. 

Theorem  3:  If  a  is  the  negative  eigenvalue  of  K  -  A(Mj  of  smallest  magnitude,  then  given 
and  two  real  spectral  shifts  A2  >  Ai,  which  are  not  eigenvalues  and  which  satisfy  one  of 
the  following 

(a)  Ai  >0 

(b)  A2  <  0 

(c)  A2  >  0  and  Ai  €  (o.O) 

Then  the  number  of  eigenvlaues  on  [Aj,A2)  is  the  difference  between  the  number  of 
positive  members  of  the  signatures  of  K  -  Ai[M|  and  K  -  A2IMI. 

Remarks: 

(1)  The  proof  of  Theorem  3  follows  directly  from  the  canonical  form,  equation  (22),  which 
has  been  shown  to  hold  for  all  stable  lumpings. 

(2)  Theorem  3  says  that  a  modified  Sturm  property  holds  for  /C  -  A[M],  but  two  shifts  are 
required  to  begin. 

(3)  The  excluded  case  is  A 2  >0,  and  Ai  <  a.  This  can  be  avoided  by  using  Ai  =  -e  >  o 
and  varying  A2  >  0  to  find  the  positive  eigenvalues.  If  the  negative  eigenvalues  are  needed, 
take  A2  =  -c  and  vary  Aj  <  0.  In  view  of  Theorem  2,  the  chances  that  a  “safe”  e  <  -a 
cannot  be  guessed  a  priori  are  negligible. 
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(4)  The  number  of  positive  members  of  the  required  signatures  can  be  computed  from  the 
number  of  sign  agreements  of  the  minors  [6j. 

5.  ALGORITHMS  FOR  THE  LUMPED  EIGENPROBLEM 

For  vibration  analysis  and  dynamic  analysis  by  modal  synthesis,  only  the  lower  portion 
of  the  spectrum  is  needed.  We  have  already  argued  that  this  portion  of  the  spectrum  is  as 
accurate  in  the  optimally  lumped  case  as  the  consistent  case.  However,  the  indefiniteness 
and/or  .singularity  of  the  mass  matrix  does  affect  the  choice  of  algorithm  for  the  solution 
of  the  eigenproblem.  Clearly  the  inverse  and  or  square  root  of  \M\  is  no  longer  available. 
Lumping  in  vibration  analysis  always  saves  the  storage  of  the  upper  triangle  of  jMj,  which 
may  be  significant,  but  to  take  full  advantage  of  lumping,  solution  algorithms  should  be 
chosen  which,  apart  from  being  indifferent  to  the  indefiniteness  and/or  singularity,  display 
a  reduced  operation  count  owing  to  the  diagonal  form. 

5.1  Algorithms 

We  have  chosen  two  algorithms  which  are  indifferent  to  the  indefiniteness  of  [M]  and 
each  of  which  saves  essentially  a  matrix-vector  multiplication  per  iteration  owing  to  the 
diagonal  form  of  the  mass  matrix.  These  algorithms  have  not  been  extensively  tested,  but 
preliminary  numerical  tests  suggest  that  they  are  as  robust  as  their  consistent  counterparts. 
In  each  case,  we  assume  that  h  <  ho  for  ho  suflficiently  small  to  guarentee  linear  divisors. 
As  pointed  out  earlier,  this  has  be  found  to  be  no  real  restriction  in  practice,  as  the  divisors 
are  usually  linear  with  few  or  even  one  element.  When  the  divisors  are  linear,  one  may 
easily  generalize  the  standard  convergence  proofs  for  each  of  the  algorithms  given  below. 

(1)  Inverse  iteration 

(|/f|-M(A<r]){IV.,.,}  =  -„[A3tl{W,)  (23) 

where  {IFfc}  is  the  current  eigenvector  iterate  and  7^  is  an  approriate  normalizing  factor 
[7].  Iteration  for  one  eigenvector  at  a  time,  as  implied  by  equation  (23)  is  not  the  usual  fi¬ 
nite  element  practice.  Instead  blcKk  iteration  is  often  employed  [6].  This  method  is  greatly 
enhanced  by  the  use  of  generalized  Jacobi  iterations  on  the  pencil  jK fc)  -  AjAffc]  obtained  by 
transformation  using  the  current  eigenvector  iterates  at  each  step.  The  transformed  pencil 
is  not  sparse,  but  approaching  closer  to  a  diagonal  one  at  each  block  iteration;  the  gener¬ 
alized  Jacobi  method  is  more  attractive  than  it  would  be  otherwise  in  this  circumstance. 
If  IMj  is  rendered  indefinite,  this  could  cause  problems  if  the  indefiniteness  carries  over  to 
the  transformed  matrix,  because  it  does  not  appear  that  Jacobi  iteration  will  work  if  [A/*] 
is  indefinite.  We  are  currently  testing  the  following  strategy;  shift  j/f)  to  make  it  and  each 
\Kk\  positive  definite,  if  necessary,  and  reverse  the  role  of  \Kk\  and  jM*]  during  the  Jacobi 
iterations.  The  reciprocals  of  the  eigenvalues  thus  obtained  are  used  where  appropriate. 

(2)  Factored,  shifted  Lanezos: 

Again,  if  neccessary,  we  shift  to  make  [K]  positive  definite  and  exchange  the  roles 
of  \K]  and  [M]-  The  Lanezos  algorithm  is  applied  in  factored  form  [8,9j  to  produce  a 
tridiagonal  matrix  with  rows  {0 . . .  7,_  i  a,  0,^  1 ...  0}: 


\A\  =  [K\^Xo[M]  \L\{L\^  =  [A] 

^.iL’{V.+  ,}=  \M  \L\~'^  m  -  a,\L]{V,}  -  3,\L]{V,}  (3,  =  0) 


(24) 


where  the  Cholesky  decomposition  is  used  to  obtain  [L].  The  normalization  factor,  7,-] 
can  be  made  equal  to  by  suitable  recursive  definition  of  a,  and  /?,,  so  that  the  resulting 
tridiagonal  matrix  is  symmetric.  The  tridiagonal  ordinary  eigenproblem  can  be  solved  any 
number  of  ways  [lO],  and  the  V,  can  be  saved  to  reconstruct  the  original  eigenvectors. 
The  attractive  feature  of  this  implementation  is  that  the  Lanczos  algorithm  produces  all 
eigenvalues  in  N  iterations  in  exact  arithmetic,  but  acts  as  an  iterative  algorithm  for  the 
extreme  eigenvalues  and  vectors  reqiiiring  very  many  fewer  than  N  iterates.  The  extreme 
eigenvalues  at  the  high  end  of  the  spectrum  for  equation  (24)  are  the  reciprocals  of  the 
smallest  positive  eigenvalues  of  the  original  shifted  pencil,  and  are  precisely  the  modes  we 
are  interested  in  getting. 

There  appear  to  be  several  other  indefinite-indifferent  algorithms  which  save  time 
if  the  mass  matrix  is  diagonal:  one  is  mentioned  in  ref.  1.  That  is  the  direct  conjugate 
gradient  minimization  of  the  square  of  the  Rayleigh  quotient  |l,]  li.  This  is  nearly  identical 
in  form  to  the  minimization  of  the  Rayleigh  quotient  itself  (the  difference  only  appears  in 
the  scaling  of  the  search  direction,  which  has  no  practical  effect).  In  fact,  a  code  devised 
for  the  consitent  case  can  be  modified  to  take  advantage  of  the  diagonal  form  of  \M\  and 
will  work  otherwise  unmodified,  even  with  zero  and/or  negative  masses  —  as  long  as  some 
care  is  exercized  in  choosing  the  initial  guess.  In  a  number  of  tests  carried  out  by  the  first 
author,  the  lumped  algorithm  was  unaffected  by  the  poles  in  the  Rayleigh  quotient,  which 
are  present  when  the  lumped  mass  matrix  is  indefinite.  The  algorithm  is  nearly  twice  as 
fast  —  half  the  number  of  matrix-vector  multiplications  are  required  per  iteration  when 
compared  to  the  consistent  mass  matrix  algorithm.  Furthermore,  the  lumped  algorithm 
often  converged  in  fewer  steps  than  the  consistent  algorithm.  The  algorithm  could  be 
defeated  only  by  choosing  a  singular  vector  or  a  vector  entirely  deficient  in  the  lower  modes 
as  an  initial  iterate.  This  situation  can  be  entirely  avoided  by  using  physical  intuition  as 
a  guide  in  the  choice  of  Fq. 

6.  TRANSIENT  ANALYSIS 

In  this  section  we  will  demonstrate  the  feasibility  and  attractiveness  of  transient  anal¬ 
ysis  by  direct  time  integration  using  optimally  lumped  mass  matrices.  Attention  will  be 
restricted  to  mass  matrix  lumpings  resulting  in  positive  and  zero  nodal  masses*  and  we 
will  consider,  as  a  model  finite  element,  the  6-node  quadratic  triangle  in  which  the  vertex 
nodes  have  zero  mass;  see  Figure  1. 

6.1  Partitioned  Mixed  Integration  Schemes 

A  time  integration  scheme  is  constructed  by  employing  a  nodal  partition  consisting  of 
implicit  and  explicit  nodal  groups;  the  implicit  nodal  group  contains  all  of  the  zero  mass 
nodes  and  the  explicit  group  contains  all  of  the  positive  m^lss  nodes.  The  partitioned, 
linear  undamped  equations  of  motion  become 


*  Research  into  the  stability  and  accuracy  of  direct  time  integration  methods  for  finite 
element  models  with  indefinite  mass  matrix  lumpings  is  currently  in  progress. 
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where  {d^}  and  {a^}  are  (n  x  1)  vectors  and  represent  implicit  approximations  to  the 
displacement  and  acceleration,  respectively,  and  {d^}  and  {a^}  are  (m  x  1)  vectors  and 
represent  explicit  approximations  to  the  displacement  and  acceleration,  respectively.  In 
equation  (25)  and  in  what  follows,  there  often  are  corresponding  statements  on  the  global 
and  element  levels,  while  some  statements  —  such  as  those  involving  global  matrix  inverses 
-  make  sense  only  on  the  global  level.  Since  the  algorithm  presented  in  this  section  performs 
as  many  operations  as  possible  at  the  element  level,  we  have  opted  to  use  the  lower  case 
symbols  indicating  element  level  expressions  whenever  a  statement  has  both  global  and 
element  level  versions,  reserving  upper  case  symbols  for  those  statements  which  make  sense 
only  on  the  global  level.  The  explicit  mass  matrix,  jm^],  is  positive  definite  and  the  stiffness 
matrix  given  by  the  assemblage  of  the  four  partitioned  submatrices  shown  in  equation  (25) 
is  assumed  to  be  positive  semi-definite  on  the  element  level,  and  its  assembled  global 
counterpart  is  assumed  to  be  positive  definite  because  of  restraints.  This  latter  restriction 
on  the  global  stiffness  matrix  does  not  seem  to  be  essential  but  simplifies  the  analysis.  The 
external  force  vector  is  partitioned  as  {p^}  and  {p^},  the  components  of  which  are  given 
on  the  element  level  by 

Pj=  [  N,<l>,(x)dn  (26) 

where  Nj  is  the  shape  function  for  node  I  of  element  e,  d>t(x)  is  either  a  body  force  or  surface 
traction  in  spatial  direction  t  and  0*  represents  element  surface  or  volume  depending  upon 
whether  is  a  traction  or  body  force,  respectively;  J  is  the  degree-of-freedom  number  of 
d.o.f.  e  at  node  /.  By  using  the  same  quadrature  rule  to  evaluate  equation  (26)  as  that 
used  in  the  mass  matrix  lumping,  all  of  the  external  loads  are  lumped  into  {p^}  and  {p^} 
is  null. 

Numerous  combinations  of  implicit  and  explicit  formulas  are  possible  for  the  time  in¬ 
tegration  of  equation  (25).  However,  because  the  implicit  mass  matrix  is  null,  the  accuracy 
and  stability  of  the  time  integration  is  directly  dependent  upon  the  explicit  method  that 
is  employed  and  is  independent  of  the  implicit  method  that  is  used  which  simply  serves  to 
define  the  average  acceleration  and  velocity  of  the  massless  nodes  whose  position  in  time 
are  exclusively  dependent  upon  the  position  of  the  adjacent  positive  mass  (explicit)  nodes. 
In  particular,  we  will  employ  the  explicit  Newmark  formula  [12j  with  /3  =  0  and  'y  =  1/2, 
which  is  similar  to  central  difference,  and  the  implicit  Newmark  formula  with  /?  =  1/4  and 
■y  —  1/2,  which  is  identical  to  the  trapezoidal  rule  (we  omit  braces  enclosing  nodal  vectors 
in  what  follows,  whenever  it  is  clear  from  the  context  to  do  so): 


<^n+l  =dn  + 


At^ 


al  +  a 


n+lJ 


*’n 


2  i«n-«n- 


(29) 

(30) 


where  subscript  n  denotes  time  nAt. 

Combining  equation  (27)  with  (25)  and  noting  that  is  null  if  the  external  forces  are 
lumped  using  the  same  quadrature  rule  as  that  used  in  the  mass  matrix  lumping,  provides 


w, . )  =■  -  i*'®i  {4  -  ^‘4  -  ^4}  pi) 

{4^,}  =  {pL,}-  i‘^'l  {<*.}  -  [*^^1  {«lf  -  Al-if  +  }  P2) 


The  computational  procedure  is  shown  in  Figure  2;  for  linear  problems,  its  implementation 
in  analysis  software  with  existing  stiffness  method  static  analysis  and  explicit  method 
transient  analysis  capabilities  is  straightforward. 

Remarks: 

(1)  Because  the  method  is  partially  implicit,  it  is  necessary  to  form  and  factor  the  global 

stiffness  matrix  However,  through  a  suitable  node  numbering  scheme,  this  matrix 

has  the  same  size  and  band  structure  as  that  associated  with  a  mesh  of  an  equal  number 
of  constant  strain  triangle  finite  elements.  Thus,  the  requirements  for  storage  and  factor¬ 
ization  are  considerably  less  than  those  usually  associated  with  quadratic  triangle  finite 
elements. 

(2)  In  the  implementation  shown  in  Figure  2,  two  internal  force  vector  evaluations  per 
time  step  are  required.  While  these  evaluations  represent  the  major  source  of  effort  per 
time  step,  a  subsequent  example  problem  will  illustrate  that  it  is  still  very  competitive  in 
terms  of  efficiency. 

(3)  These  internal  force  evaluations  can  be  carried  out  at  the  element  level,  as  pointed  out 
by  the  correspondence  implied  between  equations  (14)  and  (15). 

6.3  Stability 

An  important  consideration  regarding  the  proposed  technique  is  its  numerical  stabil¬ 
ity.  Because  an  explicit  method  is  being  employed  for  the  motion  of  the  positive  mass 
nodes,  the  procedure  will  be  conditionally  stable.  With  typical  finite  element  approxima¬ 
tions  having  positive  definite  mass  matrices,  the  stability  limit  for  explicit  computation  is 
proportional  to  the  reciprocal  of  the  highest  natural  frequency  of  the  semi-discretization 
|13|  which  is  strictly  mesh  dependent  and  has  no  physical  meaning.  In  fact,  a  substantial 
portion  of  the  upper  frequency  spectrum  is  physically  meaningless.  The  effect  of  a  mass 
quadrature  lumping  rule  that  produces  positive  and  zero  nodal  masses  is  to  force  the  mean¬ 
ingless  portion  of  the  frequency  spectrum  to  infinite;  in  the  proposed  scheme  the  motion 
at  these  frequencies  is  treated  implicitly  which  renders  their  treatment  unconditionally 
stable.  Furthermore,  because  the  remaining  portion  of  the  frequency  spectrum  has  been 
reduced,  we  expect  an  improvement  in  the  time  step  permitted  for  explicit  computation. 


14 


1.  set  initial  conditions,  n  =  0,  {£>0}  ,  {^o}  {Ao} 

2.  compute  {Dt,}  =  (D^)  -  At  {V®}  +  ^  {A^j 


3.  compute 


4.  compute 

5.  compute 


{l3'„}  =  -iA-"r'  [K'^KdL,} 

compute  element-by-element 

{At,}- 

{Pi,}  -  iK^^]  {Dl,}  -  {Dl,} 

Nib  II  im  ■■  / 


compute  element-by-element 


6.  compute  {V^,}  =  {V„^}  +  ^  [{>1^}  +  {>1^+,}] 

7.  compute  {V^j}  =  {Vf}  -H  ^  [{24^}  + 

8.  n  n  -t-  1,  go  to  2 


FIGURE  2 

Computational  procedure;  steps  4  and  6  can  be  omitted  if  the 
velocities  and  accelerations  of  the  implicit  nodes  are  not  of 
interest. 


We  define  a  stable  method  as  one  in  which  a  norm  of  the  solution,  5,  exists  at  time 
nA<  such  that 

5„  <  c  (33) 

for  all  n  where  c  is  a  nonnegative  finite  constant  jl3j;  any  vector  whose  norm  satisfies 
equation  (33)  is  said  to  been  bounded.  In  our  analysis  of  stability,  we  will  employ  an 
energy  method  similar  to  [14)  in  which  the  conditions  required  for  the  satisfaction  of 
equation  (33)  are  deduced.  The  following  definitions  are  employed 


rf-<in-i-<(.  J3^J 

d  —  <fn+ 1  dfi 

with  similar  definitions  for  the  velocity  and  acceleration  difference  and  sum. 

Moving  to  the  global  level,  the  difference  of  the  equation  of  motion  at  times  (n+  l)Af 
and  nAt  provides 

j^JE]  J  -  =  {0}  (35) 

•{•4^}  -  {i)^}  ^  =  {0}  (36) 

where  we  note  that  it  is  sufficient  to  consider  the  homogeneous  form  of  equation  (25) 
under  our  assumption  that  the  external  forces  are  bounded.  Premultiplying  equation  (35) 
by  {V^}"^  and  combining  with  the  explicit  formula,  Eqs.  (27)  and  (28)  provides 

^  iA"jl  -  {v'^}^[/i”]  {P^}  =  0  (37) 


where 


[K-\  =  [K^^]  -  [K^^]  ([K^^])  [K^^] 

Upon  rearrangement,  equation  (37)  becomes 


where 


{.4f}  \M-\{A^,]  +  {V^Y\K-\[V^) 


IM-|= 


Theorem  4:  If  \M']  and  [K']  are  positive  definite,  then  and  are  bounded. 

Proof:  Using  induction  on  equation  (39)  provides 

and  by  virtue  of  equation  (28),  the  theorem  follows.  I 

Corollary:  and  Dl^^^  are  bounded. 

Proof:  The  boundedness  of  results  from  equation  (27)  at  the  global  level,  and  then 

the  boundedness  of  results  from  equation  (31)  in  which  it  is  noted  that  a  sufficient 

condition  for  nonsingular  \K^^]  is  that  the  global  stiffness  matrix  is  nonsingular.  I 

Thus,  it  only  remains  to  determine  the  conditions  under  which  the  matrices  \K']  and 
M'  are  positive  definite. 


^  .  V  *w,  l/v  . 


Theorem  5:  \K']  =  \K^^]  -  is  positive  definite  for  arbitrary  parti¬ 

tions  if  the  assembled  stiffness  matrix 


K  = 


j^EE 


is  symmetric  and  positive  definite. 

Proof:  Pre-  and  post-multiply  equation  (42)  by  a  partitioned  vector  where 

{X^}  is  (n  •  1)  and  is  nontrivial  and  {X^}  is  (m  x  1)  and  is  arbitrary 

{X'}^  K'^]  (X'}  +  2{X^}^  [K^^]  {X®}  +  {X^}’’  [K^^]  {X^}  >  0  (43) 

For  the  particular  choice  {X^}  =  -\K^^]~^\K^^]{X^},  equation  (43)  provides 

{X®}^  [[X^^]  -  [K^^]  ([X^^])"'  [X^^]]  {X^}  >  0  (44) 

and  the  conclusions  follow.  • 

Thus,  the  stability  of  the  method  depends  exclusively  upon  the  positive  definiteness 
of  the  matrix  \M"]  =  [M^\  -  ^[X"].  If  (wf)^  and  0*,  k  =  l,2,...,n,  are  the  eigenvalues 
and  vectors  corresponding  to 

(lK-|-(wf)“|A/*]){V.»)  =  {0}  (45) 

then  [M’]  can  be  uncoupled  to  give  the  time  step  size  condition  for  stability 


where  we  refer  to  as  the  explicit  frequency,  even  though  it  is  dependent  upon  the 
implicit  part  of  the  stiffness  matrix  as  shown  in  equation  (38). 

Remarks: 

(1)  The  extreme  frequency  entering  in  the  stability  criterion,  equation  (46),  is  based  upon 
a  stiffness  matrix  that  is  the  difference  of  two  positive  definite  matrices,  as  shown  in 
equation  (38),  and  a  mass  matrix,  which  because  of  the  quadrature  lumping  rule,  has 
entries  of  larger  magnitude  than  those  associated  with  more  conventional  positive  definite 
mass  matrix  lumpings.  Thus,  the  proposed  technique  results  in  a  more  generous  time  step 
restriction  than  that  associated  with  the  exclusive  use  of  an  explicit  method  with  a  positive 
definite  mass  matrix. 

(2)  Because  the  eigenproblem  given  by  equation  (45)  involves  a  global  matrix  inversion, 
the  extreme  eigenvalue  cannot  be  bounded  by  the  maxima  of  all  unconstrained  element 
frequencies  (15].  However,  because  the  second  term  of  equation  (38)  is  positive  definite,  it 
is  possible  to  bound  the  frequency  entering  in  equation  (46)  by  the  extreme  frequency  of 
the  eigenproblem 


(^[K^^]  -  (uf  )^  [M^  ^  {ti^k}  =  {0}  (47) 

in  which  it  is  noted  that  the  eigenvalue  separation  theorem  15,16  is  applicable. 

(3)  Because  of  the  definition  of  i/i"  i,  equation  (38).  and  the  particular  nodal  partition 
that  we  are  employing  (i.e..  the  implicit  nodal  group  does  not  occupy  a  contiguous  region 
of  space,  but  rather  is  interspersed  with  the  explicit  nodal  group  such  that  each  pair  of 
implicit  nodes  is  separated  by  an  explicit  node)  the  time  step  restriction  given  by  equation 
(46)  shows  that  stability  is  not  independent  of  the  part  of  the  stiffness  that  is  treated 
implicitly.  For  example,  a  stiff  spring  placed  between  any  two  implicit  nodes  in  a  mesh  of 
quadratic  triangle  finite  elements  will  adversely  affect  the  timestep  restriction. 

6.2  Example  Problem 

In  the  examples  that  follow,  we  will  demonstrate  the  cost  effectiveness  of  the  implicit- 
explicit  optimally  lumped  quadratic  displacement  triangle  element  by  comparison  with  an 
explicit  analysis  employing  a  leading  competitor  element,  the  9-node  quadratic  displace¬ 
ment  quadrilateral  with  optimally  lumped  mass  matrices  that  are  positive  definite.  Also, 
the  loss  of  accuracy  associated  with  a  particular  ad-hoc  lumping  rule  when  applied  to 
the  quadratic  triangle  element  will  be  demonstrated  by  example.  All  computations  were 
performed  on  a  Harris  800  computer  with  double  precision  arithmetic  (48  bits  per  floating 
point  word). 

Shown  in  Figure  3  is  a  model  of  a  cantilever  beam  consisting  of  40  quadratic  displace¬ 
ment  triangular  finite  elements;  the  material  parameters  are  also  listed  in  the  figure.  The 
beam  is  unloaded  and  the  initial  condition  consists  of  the  beam  at  rest  except  for  the  free 
end  of  the  beam  (nodes  31-33,  104  and  105)  which  has  a  vertical  unit  step  velocity. 

Because  optimal  mass  lumping  is  employed,  nodes  1-33  have  zero  mass  and  are  treated 
implicitly  and  nodes  34-105  have  positive  mass  and  are  treated  explicitly.  Note  that  the 
node  numbering  scheme  is  such  that  the  size  and  band  structure  of  \K^^\  is  minimized. 
Furthermore,  because  the  explicit  matrices  \K^^\,  [K^^]  and  [K^^]  are  never  assembled, 
as  shown  in  Figure  2,  the  efficiency  of  the  explicit  pheise  of  the  computation  is  independent 
of  the  node  numbering  scheme. 

The  maximum  time  step  for  stable  computation  based  on  an  eigenanalysis  of  the  global 
matrix  equations,  Eqs.  (45)  and  (46),  is  At  <  1.93  x  10~®  sec.  Because  the  eigenanalysis 
required  to  arrive  at  this  result  is  rather  involved,  we  also  consider  a  more  conservative  and 
much  easier  to  compute  bound  based  on  an  eigenanalysis  of  the  explicit-explicit  stiffness 
matrix  only,  equation  (46),  for  a  single  element;  this  provides  At  <  1.84  x  10~®  sec.  which 
is  conservative  by  about  5%. 

Simulations  were  performed  for  a  duration  of  a  little  more  than  one  fundamental  period 
of  vibration  and  at  various  time  step  sizes  to  validate  the  stability  criteria  determined 
above;  computer  execution  times  for  these  analyses  are  listed  in  Table  1. 

Results  for  the  two  lower  time  step  sizes  considered  were  stable  and  virtually  identical; 
the  vertical  displacement  time  history  of  the  beam  tip  at  the  midsurface  (node  32)  is  shown 
in  Figure  4  for  the  simulation  employing  the  larger  step  size.  At  =  1.9  x  10“®  sec.  In  order 
to  verify  the  theoretical  stability  limit  for  the  method,  the  simulation  was  repeated  using 
A<  =  2.0  X  10"*^  sec.  and  indeed,  the  computation  was  unstable. 
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FIGURE  3 

Finite  element  models  consisting  of  40  quadratic  displacement 
triangle  elements  (Fig.  a)  and  20  biquadratic  isoparametric 
elements  (Fig.  b);  material  properties:  modulus  of  elasticity  = 
30x  10®  pst.,  Poisson’s  ratio=0.3  and  density=7.4x  10“^  . 


To  show  that  the  method  is  accurate  and  cost  effective,  comparison  simulations  were 
performed  employing  a  mesh  of  20  bi-quadratic  displacement  quadrilateral  finite  elements 
with  optimal  mass  lumping;  see  Figure  1.  Because  the  Newton-Cotes  formula  employed 
in  the  mass  quadrature  has  positive  weights,  the  resulting  mass  matrix  is  positive  definite 
and  a  strictly  explicit  method  can  be  used  for  the  time  integration.  This  is  a  particularly 
fair  comparison  because  one  quadrilateral  element  has  the  same  number  of  degrees-of- 
freedom  and  a  similar  displacement  field  as  a  macroelement  consisting  of  two  quadratic 
triangular  finite  elements.  Based  upon  an  eigenanalysis  of  a  single  element,  the  time  step 
restriction  was  determined  to  be  At  '  1.65  »  10  ®  sec.  Simulations  were  performed  using 
various  time  step  sizes  and  for  a  duration  equivalent  to  that  employed  for  the  triangular 
finite  element  computations:  the  computer  execution  times  are  given  in  Table  1  and  the 
vertical  displacement  time  history  of  the  beam  tip  at  the  mid-surface  (node  32)  for  the 
largest  stable  time  step.  At  -  1.7  >  10  *'  sec.,  is  shown  in  Figure  4  beside  the  results 
for  the  implicit-explicit  optimally  lumped  triangle  element.  These  results  agree  extremely 
well  except  for  a  slight  period  error.  Furthermore,  a  comparison  of  the  execution  times 
for  the  simulations  surprisingly  show  the  implicit-explicit  triangular  element  analysis  to  be 
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Quadratic  Triangle  with  Optimal  Mass  Lumping 


At{  X 10  ®  sec. ) 

#  steps 

execution  time  (minutes) 

1.8 

1775 

94 

1.9 

1685 

89 

2.0 

1600 

unstable 

Quadratic  Quadrilateral  with  Optimal  Mass  Lumping 

At  (  <  10  ^  sec. ) 

#steps 

execution  time  (minutes) 

1.6 

2000 

137 

1.7 

1800 

129 

1.8 

1 775 

unstable 

Quadratic  Triangle  with  Ad-hoc  Mass  Lumping 

At{  X  10~^'  .see. )  Esteps  execution  time  (minutes) 
1.2  2665  76 

TABLE  1 

Time  step  sizes  and  computer  execution  times  for  example 
problem. 


substantially  more  economical  than  the  explicit  quadrilateral  element  analysis.  The  reason 
for  this  efficiency  is  threefold:  first,  the  effort  required  to  store  and  factor  \K^^]  was  minimal 
for  this  problem.  In  fact,  the  size  and  band  structure  of  this  matrix  is  identical  to  that  for  a 
mesh  of  40  constant  strain  triangular  elements.  Although  the  number  of  equations  in  \K‘^] 
is  considerably  smaller  than  the  total  number  of  equations  (in  this  example,  it  h8is  about 
30%  of  the  total  number  of  equations)  it  is  noted  that  its  storage  and  factorization  become 
increasing  more  inconvenient  with  increasing  problem  size.  The  second  source  is  the  fact 
that  fewer  time  steps  are  required  because  of  the  improved  stability,  thus,  approximately 
6%  fewer  steps  were  needed.  The  last  source  of  savings  stems  from  the  cost  associated  with 
internal  force  evaluations.  Even  though  the  implicit-explicit  triangle  algorithm  requires 
two  internal  force  evaluations  per  time  step;  see  Figure  2,  as  opposed  to  one  for  a  strictly 
explicit  analysis,  the  fact  that  only  six  integration  points  per  two-triangle  macro  element 
are  required  compared  to  nine  integration  points  for  the  quadrilateral,  coupled  with  the 
fact  that  the  triangle  evaluations  are  easier  (fewer  multiplications)  renders  the  triangle 
internal  force  evaluations  slightly  more  efficient  per  time  step. 

In  order  to  demonstrate  the  importance  of  using  an  integration  rule  of  sufficient  ac¬ 
curacy  in  the  mass  quadrature,  we  will  compare  the  results  obtained  by  the  previous 
implicit-explicit  optimally  lumped  triangle  simulation  with  an  explicit  analysis  employing 
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FIGURE  4 

Cantilever  beam  tip  displacement  time  histories  for  quadratic 
triangle  and  biquadratic  isoparametric  element  analyses  em¬ 
ploying  optimal  mass  lumping. 


an  identical  mesh  of  quadratic  triangle  elements  with  ad-hoc  mass  lumping  to  preserve 
the  positive  definiteness  of  the  mass  matrix.  The  particular  scheme  that  was  employed  is 
detailed  in  [3,  and  consists  of  letting  the  diagonal  entries  of  the  lumped  matrix  equal  the 
diagonal  entries  of  the  consistent  mass  matrix  multiplied  by  the  ratio  m/s  where  m  is  the 
total  mass  of  the  element  and  s  is  the  sum  of  the  diagonal  entries  of  the  consistent  mass 
associated  with  translational  degrees-of-freedom  in  one  coordinate  direction  only.  The  re¬ 
quired  computer  execution  time  is  given  in  Table  1  and  the  vertical  displacement  time 
history  of  the  beam  tip  at  the  mid-surface  is  shown  in  Figure  5.  Although  the  phase  of 
the  two  solutions  agree  quite  well,  the  amplitudes  disagree  by  as  much  as  35%.  Therefore 
it  is  quite  apparent  that,  at  least  for  the  particular  ad-hoc  lumping  rule  employed  in  this 
simulation,  that  substantial  inaccuracies  can  result  from  non-optimal  mass  lumping  rules. 

7.  CONCLUSIONS 

We  have  demonstrated  that  optimal  lumping  of  masses  for  higher-order  elements  can 
be  an  viable  computational  technique,  in  second-order  problems.  When  only  zero  masses 
result  from  such  lumping,  we  have  shown  that  both  vibration  and  dynamic  analyses  can 
be  carried  out  in  a  manner  which  is  more  accurate  than  an  ad-hoc  lumping  which  pre- 
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FIGURE  5 

Cantilever  beam  tip  displacement  time  histories  for  quadratic 
triangle  element  analyses  employing  optimal  mass  lumping  and 
ad-hoc  mass  lumping. 


serves  positive  mass,  and  even  more  computationally  efficient  than  a  widely-used  scheme 
of  comparable  accuracy  which  has  positive  masses.  Implicit  treatment  of  the  resulting 
infinite  frequencies  appears  to  be  a  favorable  trade-off  in  dynamic  analysis,  as  long  as  the 
number  of  zero  masses  is  not  excessive.  We  envisage  that  most  of  the  application  of  the 
techniques  decribed  here  will  be  in  structural  and  elzisto-dynamic  problems,  but  it  may 
be  of  interest  to  mention  that  the  motivating  example  in  devising  the  dynamic  algorithm 
presented  here  is  that  of  fluid  mechanics.  There  the  pressure  is  in  somes  sense  analogous 
to  the  modes  of  infinite  frequency  introduced  by  zero  masses.  It  has  long  been  known  that 
the  pressure  must  be  treated  implicitly  in  algorithms  which  treat  the  velocity  explicitly. 
The  reader  may  have  observed  that  there  is  a  subtle  difference  with  the  pressure  in  incom¬ 
pressible  fluid  dynamics,  however.  The  pressure  at  the  current  iterate  of  a  velocity-explicit 
Navier-Stokes  method  must  be  chosen  to  make  the  velocity  field  weakly  divergence-free 
at  the  next  iterate.  The  situation  here  is  simpler;  the  massless  nodes,  involving  modes 
of  infinite  frequency,  are  carried  by  the  stiffness  of  connected  nodes  with  non-zero  mass. 
As  a  consequence,  the  motion  of  the  nodes  with  nonzero  mass  is  computed  first,  and  the 
motion  of  the  massless  nodes  is  computed  from  the  motion  of  the  nonzero  masses.  Such  a 
computational  ordering  would  by  unconditionally  unstable  in  fluid  mechanics. 
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In  vibration  analysis,  the  infinite  frequencies  introduced  by  lumping  can  be  viewed 
as  the  setting  of  some  of  the  inaccurate,  high  frequency  portion  of  the  discrete  spectrum 
to  infinite  frequency.  Heuristically,  this  is  similar  to  what  is  done  in  static  condensation, 
which  is  equivalent  to  creating  infinite  frequencies  an  a  more-or-less  ad-hoc  basis.  Static 
condensation  is  likely  to  lead  to  reduced  accuracy,  whereas  optimal  lumping  is  guided  by 
'h'orons  notions  of  stable  and  accurate  quadrature,  and  accuracy  will  not  be  lost.  Of 
(oiiise.  optimal  hm.ping  will  not  accomplish  the  same  goal  as  static  condensation,  which 
is  to  replace  all  but  a  few  of  the  most  accurate  vibration  modes  by  ones  with  infinite 
Irequency.  thereby  greatly  reducing  the  dimension  of  the  eigenproblem.  It  nevertheless 
appears  that  optimal  lumping  can  be  used  in  a  similar  spirit  to  reduce  the  computational 
cost  of  eigenanalysis.  Optimal  lumping  may  also  be  viewed  as  a  means  of  pushing  to  the 
limit  the  notion  of  modifying  the  upper  portion  of  the  discrete  spectrum  without  suffering 
an  accuracy  loss. 

This  finally  brings  us  to  the  question  of  negative  masses.  In  spite  of  first  appearences 
to  the  contrary,  we  have  shown  that  they  do  not  interfere  with  vibration  analysis.  We 
like  to  think  of  the  negative  masses  as  ones  that  have  been  forced  to  even  “higher  than 
infinite”  frequency  by  the  lumping  scheme,  though  this  may  only  be  a  somewhat  quaint 
descriptive  analogy.  There  should  be  a  very  real  doubt  about  the  dynamical  viability  of 
negative  masses:  they  introduce  high  frequency,  unconditionally  unstable  modes  into  the 
problem,  which  seems  the  last  thing  anyone  would  want  to  do  in  dynamic  analysis.  We  do 
not  believe  it  is  as  bad  as  that;  the  saving  grace  may  lie  in  the  fact  we  have  established 
here:  the  negative  masses  give  rise  to  very  high  frequency  modes.  This  portion  of  the 
spectrum  is  just  the  one  subject  to  algorithmic  damping  by  many  good  integrators.  To 
put  it  another  way,  many  ODE  integrators  have  stable  regions  which  include  substantial 
portions  of  the  right  half  plane  (backwards  Euler’s  method,  for  one,  though  we  are  not 
recommending  it  as  a  practical  method).  We  are  currently  invesigating  the  possiblity  of 
devising  time-stepping  procedures  with  sufficient  algorithmic  damping  to  be  indifferent  to 
spurious  unstable  modes  introduced  by  optimal  lumpings  with  negative  masses. 
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APPENDIX 

Sharpening  of  G.  Fix’s  result  [oj*  on  quadrature  errors  when  a(-,  •)  =  a(",  •).  Following 
ref.  5,  element-by-element  quadrature  errors  are  of  the  fornn  (applied  to  scalar  component 
of  u**.  call  it  u^): 

//*'■  I  D  ■  (wS-M'  dV:  ’'yi  -  (A.l) 

where  '•  -  typical  derivative  multi-index,  ii,  -  degree  of  lowest  order  polynomial  not 
integrated  exactly.  We  assume,  as  does  ref.  5,  that  the  inverse  inequality 
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on  each  element.  There  are  two  cases.  (.\)  Stability:  and  are  arbitrary  trial  functions 

and  (B)  accuracy;  is  the  f.e.m.  solution  and  is  an  arbitrary  trial  function.  In  case 
(B)  for  0  <  s  <  p 

M.3) 

€ 

where  u  is  the  exact  solution.  The  left-hand  side  sums  to  ||  u*  ||^  for  s  =  0  or  1. 


Case  A:  (A.l)  is  the  sum  of  terms  of  the  form 

Ee  =  Ch^'  ^ \{D^u^)  (£>^'*)|  dV;  |a|  -h  |/?|  =  h|  =  £o 
Suppose  |a|  0  and  |/?|  ^  0,  then 
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where  |aol  =  \0o\  =  1-  Thus 

E.  <  i|u‘  !|,  ,  II II.  .  =  C.h^  II  u*  11.,,  11. -  II.,, 

Now 

(u^v^)  dV 

<  E  s  (e  II  “*  ii;,.]  ’  (e  II  Ii  II.  i!  II . 


0,e 


(A.4) 


(A.5) 


(A.6) 


*  Reference  numbers  in  the  Appendix  refer  to  the  Reference  section  of  the  main  body 
of  the  paper. 
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If  either  jai  =  0  or  jpl  =  0,  then  the  following  case  is  representative: 


A*"  /  (D-'u*)  v*  <  i;  D-Iu*  i; . , ,  ;i  V*  II , , 


(A.7) 
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and  summing  as  before  shows  in  general  that 


|6{u\v'‘) -S(u\t;'‘)|  <  Cft"* 
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where  6(-,")  refers  to  the  L2  inner  product  in  one  vector  component  of  the  trial  spaces 
defined  in  the  main  body  of  the  paper.  Since  the  same  interpolations  and  quadratures  are 
used  in  each  component,  the  result  of  equation  (A.8)  carries  over  to  6(u**,v*‘);  this  proves 
stability. 

Remarks: 

1.  Fix’s  result  can  be  sharpened  because  &(•,  •)  involves  no  derivatives.  One  derivative  “on 
each  side”  must  be  left  behind  to  compare  to  ||  •  Hj,  hence  two  fewer  powers  of  h  are  lost 
in  the  process. 

2.  When  all  £0  derivatives  are  on  one  side,  one  must  be  “moved  over”.  The  negative  norm 
argument  is  just  a  compact  way  of  applying  the  divergence  theorem  and  arguing  away  the 
boundary  integral. 

3.  The  above  argument  shows  that  any  quadrature  formula  with  degree  of  precision  >  0 
is  stable  for  6(-,  •).  Accuracy  is  another  matter. 


Case  B:  We  argue  that  (A.l)  is  small  indirectly,  making  use  of  (A. 3).  Assume  that  a 
non-conforming  triangle  of  (complete)  degree  0  <  9  <  p  can  be  inscribed  in  every  element 


of  the  mesh,  and  that  the  conditions  on  the  mesh  are  such  that  the  standard  error  analysis 
jl7j  applies  to  elements,  so  that  if  is  a  standard  approximation  to  from  the  implied 
space  of  complete  (discontinuous)  polynomials 
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Note  that  q  -  0  is  allowed  giving  only  boundedness  of  the  left-hand  side.  Further  let  us 
assume  that  b(u^.v^)  =  b{v^‘.v^)  (whereas  h{t/^‘.v^)  may  differ  from  b{u^,v^)).  Then 

b(u\v^)  -  b(u\v^)  -  6(u''  -  -  6(u''  -  u\v^)  = 

„  /  -  i'*)  </!•  -  .p  ^  £  „-j  (p*  -  vH 

f  k 

<  l!u‘  1: f*  li, ,  =  l|p‘  II V*  II, , 

€  e 

where  the  inequality  is  determined  from  the  assumption  of  stability.  The  subscript  k  refers 
to  function  value  at  the  quadrature  point.  Note  we  assume  p  =  constant;  the  result 
will  apply  when  p  is  a  smoothly  varying  coefficient  [17].  We  accumulate  the  error  as  in 
equation  (A. 6),  and  finally  arrive  at  the  desired  result: 


|6(«\u^)  -  6(u\u'*)i  <  li  u  lip^,  II  v'*  II, 


(A.IO) 


Remarks: 

1.  For  the  bilinear  element,  9  =  0,  leads  to  b(u^yV^)  =  b(u^,v^),  thus  the  entire  solution 
is  0(h)  (optimally)  accurate  [5|.  For  the  biquadratic  element,  9  =  0  or  1  leads  to  0{h^), 
which  is  also  optimal. 

2.  Bi  -  p  polynomials  have  a  power  of  h  to  spare.  With  linear,  quadratic  and  cubic 
triangles  the  exact  integration  condition  cannot  be  met  unless  9  =  0, 0  and  1  resp. 

3.  All  elements  discussed  in  the  body  of  the  paper  lead  to  recovery  of  full  accuracy  by 
allowing  p  -  9  <  2. 
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